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g ■ ABSTRACT 

All the analyses of Cosmic Microwave Background (CMB) temperature maps 
O \ up-to-date show that CMB anisotropies follow a Gaussian distribution. On the 

\ other hand, astrophysical foregrounds which hamper the detection of the CMB 

\ angular power spectrum, are not Gaussian distributed on the sky. Therefore, 

.!_H \ they should give a sizeable contribution to the CMB bispectrum. In fact, the first 

^ \ year data of the Wilkinson Microwave Anisotropy Probe (WMAP) mission have 

\ allowed the first detection of the extragalactic source contribution to the CMB 

bispectrum at 41 GHz and, at the same time, much tighter limits than before 
to non-Gaussian primordial fluctuations. In view of the above and for achieving 
higher precision in current and future CMB measurements of non-Gaussianity, 
in this paper we discuss a comprehensive assessment of the bispectrum due to 
either uncorrelated or clustered extragalactic point sources in the whole frequency 
interval around the CMB intensity peak. 

Our calculations, based on current cosmological evolution models for sources, 
show that the reduced angular bispectrum due to point sources, hps, should be 
detectable in all WMAP and Planck frequency channels. We also find agreement 
with the results on hps at 41 GHz coming from the analysis of the first year WMAP 
data . Moreover, by comparing hps with the primordial reduced CMB bispectrum, 
we find that only the peak value of the primordial bispectrum (which appears at 
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I ~ 200) results greater than bps in a frequency window around the intensity peak 
of the CMB. The amphtude of this window basically depends on the capability 
of the source detection algorithms (i.e., on the achievable flux detection limit, 
Slim, for sources). Finally, our current results show that, at low frequencies (i.e., 
u < 100 GHz), the angular bispectrum of a clustered distribution of sources 
appears not substantially different from that of Poisson distributed ones, by 
using realistic angular correlation functions suitable to be applied to the relevant 
source populations. On the other hand, we also find that at higher frequencies 
(i.e., u > 300 GHz), the clustering term can greatly enhance the normalization 
of bps. 

Subject headings: extragalactic point sources: number counts, angular distribu- 
tion - cosmic microwave background: bispectrum 



1. INTRODUCTION 

The analysis of the first year Wilkinson Microwave Anisotropy probe (WMAP) data 
clearly show that temperature fluctuations of the Cosmic Microwave Background (CMB) 
follow a Gaussian distribution (Komatsu et al. 2003; Spergel et al. 2003) in agreement with 
the standard inflation paradigm. Therefore, the third moment of this distribution and its 
angular bispectrum, Ci^i^i^, the harmonic transform of the three-point correlation function, 
should be zero and the angular power spectrum, C^, specify all the statistical properties 
of CMB anisotropy. Anyway, significant non-Gaussianity can be introduced, for instance, 
by second order relativistic effects or by features in the scalar field potential, which can 
produce an angular bispectrum potentially detectable in all-sky maps as those provided by 
the WMAP and Planck satellite missions (Bennett et al. 2003a; Mandolesi et al. 1998; Puget 
et al. 1998). 

In the standard infiationary scenario, the quantum fiuctuations of the scalar field are 
the origin of the matter and radiation fiuctuations in the Universe. In slow-roll infiation, 
weak non-Gaussian fiuctuations are generated by non-linearity in inflation (Salopek and 
Bond 1990) or by features appearing in the inflation potential (Gangui et al. 1994) which 
can produce small deviations from gaussianity in the CMB temperature fluctuations. Second 
order relativistic effects can give also rise to a detectable non-Gaussian contribution (Pyne 
and Carroll 1996). All these effects can be summarized by the following expression for the 
curvature perturbations. 
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^x)^M^)+Umx)-mx))). (1) 

where is the Gaussian part of the perturbations, /„; is the non-hnear couphng param- 

eter, given by a certain combination of the slope and the curvature of the inflaton potential, 
and means the statistical ensemble average (Falk, Rangarajan, and Srednicki 1993; Gangui 
et al. 1994; Wang and Kamionkowski 2000). In general, high values of fni can be obtained 
just with a small tilt of the spectrum and these values could give rise to a bispectrum signal 
potentially detectable in current as well as future all sky CMB maps. Moreover, the value 
of fni could be also increased by deviations from standard inflationary models (see, e.g., 
Kofman et al. 1991). 

The great astrophysical interest in the detection of primordial non-Gaussianity has 
stimulated the development of many different methods suitable in analyzing CMB anisotropy: 
Minkowski functional on the sphere (Schmalzing and Gorski 1998), statistics of excursion 
sets (Barreiro et al. 2001), wavelets (Pando et al. 1998; Barrciro ct al. 2000; Cayon et 
al. 2001), the bispectrum (Luo 1994; Komatsu and Spergel, 2001; Komatsu et al. 2002; 
Santos et al. 2002). The main result of all these analyses is that no positive detection of 
non-Gaussianity in CMB maps has been confirmed up to now: COBE DMR (Bromley and 
Tegmark 1999; Cayon et al. 2003; Komatsu et al. 2002) Maxima-1 (Wu et al. 2001; Santos et 
al. 2002) and Boomerang (Polenta et al. 2001) data appear all compatible with the Gaussian 
hypothesis. On the other hand, Ferreira et al. (1998) have measured 9 equilateral modes 

1 /2 

{ii = £2 = is) of the normalized bispectrum, Bi^i^i^/C^ , on the COBE DMR map, claiming 
detection at ii = £2 = is = 16 whereas Bromley and Tegmark (1999) find that few individual 
pixels - which could be contaminated by both foreground emission and instrumental noise 
- are responsible of most of the claimed primordial signal. 

A great challenge to precise measurements of non-Gaussianity in the CMB is set by 
astrophysical foregrounds which constitute an unavoidable limitation : see, e.g., "Microwave 
Foregrounds", 1999, A. de Oliveira Costa and M. Tegmark eds., ASP Conf. Ser. Vol.181, 
for a thorough discussion on the subject. At present, even applying the most efficient me- 
thods for component separation, the residual astrophysical signal - which is not, in general, 
Gaussian distributed - can contaminate a non negligible number of pixels giving rise to 
a detectable contribution to non-Gaussianity. Other studies have discussed the effect of 
Galactic foregrounds on non-Gaussianity (see, e.g., Komatsu et al., 2002). In this paper, 
we focus on extragalactic point sources whose contribution to CMB anisotropics has been 
analyzed in detail (Toffolatti et al. 1998; Sokasian et al., 2001; Park, Park & Ratra 2002) 
whereas their contribution to the CMB bispectrum has been poorly studied up to now. A 
first estimate of the reduced bispectrum, bps, due to Poisson distributed point sources at 
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two frequencies only (90 and 217 GHz) has been recently presented by Komatsu & Spergel 
(2001), hereafter KSOl. This paper has shown that undetected point sources, i.e. sources at 
fluxes below the detection threshold of the experiment, give rise to a non zero bps value which 
could hamper the detection of the primordial signal. Very recently, by the analysis of the 
first year WMAP observations, Komatsu et al. (2003) have claimed the first detection of the 
bispectrum due to extragalactic sources. In view of the above, it is clearly of great interest 
to perform a thorough analysis of the contribution of point sources to the CMB bispectrum. 
Therefore, we extend the calculation of KSOl to the whole frequency range of the WMAP 
and Planck missions, using model number counts able to reproduce the data available so far. 
Moreover, we take also into account the effect of correlated positions of point sources in the 
sky for studying how much source clustering can enhance the value of the reduced angular 
bispectrum, bpg. 

The outline of the paper is as follows. In Section 2 we briefly summarize the formalism 
used to estimate the reduced bispectrum and show detailed calculations of its value due 
to Poisson distributed point extragalactic sources. In Section 3, we model the clustering 
properties of sources at microwave frequencies by 2D simulations and we present the first 
estimates to date of the CMB bispectrum produced by clustered sources. Finally, in Section 
4, we discuss our main results and present the conclusions. We assume a standard Cold 
Dark Matter (CDM) cosmology (^2^=1-0, Hq — 50 km/s/Mpc) but we discuss our findings 
also in the case of a flat ACDM model with Qm — 0.3 and 0,a — 0.7. Anyway, as discussed 
below, the conclusions are largely independent of the adopted cosmological model. 

2. THE BISPECTRUM DUE TO UNCLUSTERED POINT SOURCES 

2.1. General formalism 
CMB temperature fluctuations AT(n)/T are usually expanded into spherical harmonics 



where fl denotes the full sky (in the case of an incomplete sky coverage, the integral is done 
over the observed sky area, flobs] see, e.g., Komatsu et al., 2002). The angular third moment 
of CMB temperature fluctuations is deflned as 



and the averaging is over the ensemble of realizations. If the fluctuations are gaussian 
distributed, the third moment and all the higher order moments are zero. 




(2) 




(3) 
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Due to the rotational invariance of the universe 

^™=(i\ (4) 

where Ci-^i^i.^ is the angular bispectrum and the matrix is the Wigner-3j symbol. Since we only 
have access to one sky we need an estimator of the angular bispectrum: the best unbiased 
estimator (Gangui and Martin 2000) is the angle-averaged bispectrum 

Another related quantity is the reduced bispectrum bi^^i^i^, which can be expressed easily 
in terms of the angular bispectrum 



(2/i + 1) (2/2+1) (2/3 + 1) f h h h 



Chhh = V ] X Khh (6) 

V 47r V "^1 "^2 / 

and which contains all the physical information in B^'^^^^'^^ . The reduced bispectrum can 
also be related to the third moment 

(7) 

where is the Gaunt integral 

SZ7r' = / d'{n)Yi,„,, (n) 1^,^, (n) Y^,^, (n) (8) 

The reduced angular bispectrum is a very convenient quantity, since it can be easily calcu- 
lated in the flat two-dimensional approximation (Santos et al. 2001). 

Having summarized the general formalism, we remind here the calculation of the reduced 
bispectrum expected in the case of unclustered point sources, i.e extragalactic point-like 
sources which follow a Poisson distribution in the sky. Toffolatti et al. (1998) found that 
this assumption is fairly good at microwave frequencies, if sources are not subtracted down 
to very faint flux limits, thus greatly decreasing Poisson fluctuations whereas leaving almost 
unaffected the clustering term. Under this assumption, the sky distribution of extragalactic 
sources produces a white noise spectrum and, from Eq.(7), we define 

bps^btZr = {iT-{T)r)/T' (9) 

i.e., equal to the skewness of the sky temperature distribution. Therefore, 6ps=constant at 
all scales i. 
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The reduced bispectrum, bps, can be then calculated as follows 



-'ps 

rSum J„ 

hiM,^9\x)l dSS'— (10) 







djTt 

where — is the differential source count per unit solid angle, S the flux and a the conversion 
db 

factor from fluxes to temperatures 

where x = hv/ksT. The integral in Eq.(lO) has to be computed up to the flux detection 
limit foreseen for the experiment, Sum-i since only undetected sources are contributing to the 
estimated b^g- 



2.2. Number counts of extragalactic sources 

To perform the integral in Eq. (10) we have used the differential counts corresponding to 
the cosmological evolution models for radio and far-IR selected sources discussed in Toffolatti 
et al. (1998) (hereafter T098). This choice represent a clearly better approximation than 
simple power law counts, as assumed by KSOl with the purpose of obtaining a first estimate 
of bps- However, since they used the model counts of T098 for the extrapolation of source 
counts at 94 GHz, with the flux detection threshold of Sum ~ 2 Jy (the estimated Scr 
detection hmit of WMAP W-band) , it is not surprising that their estimates result in good 
agreement with our current ones. 

As for data on radio counts, preliminary measurements at 15.2 GHz by Taylor et al. 
(2001) suggested that the extrapolation of the model counts of T098 - which successfully 
account for the observed source counts down to 5" ~ 0.1 mJy in flux and up to 8.44 GHz - 
overestimates their data by a factor of ~ 1.5 2 at the survey hmit (~ 20 mJy), implying 
that the simple assumptions on source spectra and/or on cosmological evolution could not 
have been the most correct ones. On the other hand, the recently published WMAP analysis 
of the foreground emission (Bennett et al. 2003b) has provided a full sky catalogue of 208 
bright extragalactic sources with fluxes S > 0.9 — 1.0 Jy, of which only five objects could 
be spurious identifications. ^ The whole sample give an average flat {a = 0.0 ± 0.2) energy 



^Only very few sources are detected by WMAP at fluxes S < 0.9 Jy and only two at galactic latitude 
\b"\ < 10° according to Table 5 of Bennett et al. (2003b). Moreover, at fluxes 5 < 1.2 ^ 1.3 Jy the sample 
appears to be statistically incomplete. 
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spectrum in full agreement with the assumptions of T098 and number counts of bright 
sources at 41 GHz which appear to fall below the prediction of T098 at 30 GHz by a factor 
of ~ 0.66 (but only in the faintest flux bins). Direct calculations of the 33 GHz counts by the 
source catalogue in Table 5 of Bennett ct al. (2003b) give 155 sources at Sum — 1-25 Jy on a 
sky area of 10.38 sr (|6^^| > 10°) where the sample appears to be statistically complete. By 
this sample we could estimate again the WMAP differential counts finding an average offset 
of ~ 0.75 with the T098 model predictions (see Toffolatti et al. (2003) for more details). 

Moreover, two other recently published independent samples of cxtragalactic sources at 
31 and 34 GHz - from CBl (Mason ct al. 2003) and VSA (Taylor et al. 2003) experiments, 
respectively - show that the T098 model correctly predicts number counts down to, at 
least, 5" ~ 10 mJy. Another independent check of the assumptions of the T098 model is 
coming from the few detections of extragalactic point sources in the Boomerang maps. These 
sources show flat or decreasing spectra up to 240 GHz (Coble et al. 2003) which is again in 
agreement with the population mixture of T098. Therefore, given that almost all up-to- 
date observations are confirming the predictions on source counts of the T098 model up to 
30-34 GHz, we can be very confident in applying this model for our current predictions, at 
least up to frequencies < 100 -^ 150 GHz, where "flat" -spectrum sources dominate the bright 
counts. 

On the other hand, a lot of new data on far-IR/sub-mm source counts have piled up 
since 1998, in particular by means of the SCUBA and MAMBO surveys. These current 
data arc better explained by new physical evolution models of far-lR selected sources, as the 
ones presented by Guiderdoni et al. (1998) and by Granato et al. (2001). Therefore, we also 
compare our first estimate of bps at 545 GHz with the ones obtained by directly integrating 
the differential counts of these two latter predictions. Anyway, since the reduced bispectrum 
due to undetected sources is produced by all sources below a given detection limit, the 
integral in eq. (10) results not much affected by the adopted evolution model (see below). 

In Table 1 we present our current estimates of the reduced bispectrum, bps, due to 
Poisson distributed extragalactic sources at various frequencies in the range 15-850 GHz and 
applying different fiux limits, Sum, for source detection. The frequencies and flux detection 
limits has been chosen to directly compare with those of the MAP and Planck satellite 
missions. All our estimates have been calculated by taking into account both the radio and 
far-lR source populations at each frequency channel. Anyway, at frequencies below ~ 150 
GHz the contribution of far-IR selected sources is negligible whereas at u >500 GHz the 
same applies to radio selected ones. The values of bps are, obviously, at a minimum close to 
the CMB intensity peak where bright sources reach a minimum number. 

It is worthy to remind that the flrst estimate of KSOl gave the value bps ^ 2 x 10~^^ at 
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94 GHz, with Sum ~ 2 Jy whereas our current estimate is ~ 3.7 x 10~^^; at 217 GHz and 
with Slim ~ 0.2 Jy their value was ~ 5 x 10~^^ whereas we obtain ~ 3.8 x 10~^^. The general 
agreement and, at the same time, the little mismatch with our estimates are not surprising 
since KSOl used the T098 model for estimating bps but they extrapolated T098 number 
counts with a constant slope down to faint fluxes. 

It is also interesting to compare 6ps estimates with the reduced CMB bispectrum ob- 
tained in non-Gaussian models with fluctuations of the inflation potential given by eq. (1). 
This quantity has also been carefully calculated by KSOl. They have shown that the equilat- 
eral bispectrum bm - i.e., the reduced angular bispectrum bi-^ij.^ when h = h = h ^ multiplied 
by exhibits a peak at multipole ^ ~ 200. The height of this peak is proportional to the 
factor /„;, the coupling parameter. 

In Figure 1 we also compare our estimates of bps with the primordial CMB equilateral 
bispectrum, bm, at its peak value ~ 200) by adopting fni = 100, 10, 1. As also indicated 
in Table 1, we explore the whole frequency interval around the intensity peak of the CMB. 
Looking at Figure 1, the peak value of the primordial bm appears greater than bps in a 
frequency window around the CMB intesity peak: obviously, the window sets larger by 
increasing the value of fni and by reducing the source detection limit, Sum- If fni — 100 
and with Sum — 0.1 Jy, the window spans from ~ 30 to ~ 500 GHz. On the other hand, if 
fni = 10 and Slim — 2 Jy, the window shrinks to 100 < ly < 350 GHz. The line corresponding 
to fni = 1 is plotted only as a reference, since the ideal experiment requires fni > 3, in order 
to obtain S/N > 1 (Komatsu and Spergel 2001). 

As an example, in Figure 2 we show our estimates of bps at three frequencies, 30, 44 and 
70 GHz, and for S'/jm=l,0.1 Jy comparing them with the primordial equilateral bispectrum, 
bill, calculated by assuming fni = 10. It is apparent that at this three frequencies, and, thus, 
for all frequencies below ~ 70 ^ 80 GHz, bps has always a higher value than the primordial 
CMB bispectrum if Sum > 1 Jy- On the other hand, if we are able to detect and remove 
sources brighter than Sum — 0.1 Jy, which could be attaineable by Planck (Vielva et al., 
2003), the peak value of the primordial CMB bispectrum should be detectable at frequencies 
> 40 50 GHz, for this choice of fni- We checked that fixing fni — 100, the hmits for the 
detectability of the primordial bispectrum move down in frequency to ~ 50 and ~ 20 GHz, 
if Slim = 1, 0.1 Jy, respectively. We have limited our comparisons to these values of since 
WMAP data (Komatsu et al. 2003) already imply that fni < 134, at the 95% confidence 
level. 
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Table 1: Sample values of bps in the 15-850 GHz frequency interval and for different flux 
detection limits, Sum- All the Planck channels are indicated. 
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"Central frequencies of the K-, V-, and W- MAP channels. The two remaining channels of the MAP mission, 
Ka- and Q-, have not been indicated because their central frequencies, 33 and 41 GHz, are very close to the 
ones of the Planck LFI instrument at 30 and 44 GHz. 
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Fig. 1. — Reduced equilateral bispectrum, bps, due to Poisson distributed point sources in the 
15-900 GHz frequency interval and for several flux detection limits, Sum, for sources. The peak 
value of the primordial CMB bispectrum (at i 200) generated by a quadratic potential is also 
plotted for different values of the coupling parameter, /„;. 
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As shown by KSOl, the contribution of point sources and the one coming from the 
Sunyaev-Zeldovich effect can be separated from the primordial CMB angular bispectrum via 
chi-square minimization, by taking into account the different shape of the three contributions 
and thanks to the sensitivity of satellite observations. According to Table 3 in their paper, 
the signal-to-noise ratios for bps are 2.2 x lO^^feps for WMAP and 52 x lO'^'^bps for Planck, 
respectively. The source contribution to the CMB bispectrum should therefore be detectable 
by WMAP if bps > 2 x 10~^^ (4(7 limit) and, more easily, by Planck since it should be 
sufficient that bps > 10~^^ (4c), in this latter case. A comparison of the above detection 
hmit for Planck with our current estimates of Table 1 shows that bps should be detectable 
in all Planck channels even lowering the source detection limit down to Sum ~ 0.1 Jy, 
which currently appears as an achievable limit, at least for all frequencies below ~ 200 GHz 
(Vielva et al. 2003). If we take the WMAP limit, we find that bps could be detected also at 
all WMAP frequencies (only at 1^20" level in the 94 GHz channel) if only sources not fainter 
than Slim ~ 1-0 Jy are detected and removed. This appears as a realistic possibility given 
that only very few sources in Table 5 of Bennett et al. (2003b) show fluxes S < 0.9 -r- 1.0 Jy. 

Very recently, Komatsu et al. (2003) published the first detection of bps at microwave 
frequencies: they report the value of (9.5±4.4) x 10^'^' fiK'^'si^, by adopting Sum — 0.75 Jy for 
source detection in the 41 GHz WMAP channel. Our current findings are in good agreement 
with the claim of Komatsu et al. (2003): in fact, by directly integrating the differential 
source counts of T098 up to Sum = 0.75 Jy at 41 GHz, we find bps ~ 13-e- 14 x 10~^ //X^sr^, 
where the scatter is due to the different cosmology (Qm = 1-0 or Qm — 0.3). These results 
are compatible with the detection of Komatsu et al. (2003) at la level whereas we have to 
apply a correction factor of ~ 0.7 for obtaining the best fit value of Komatsu et al. (2003). 
If we use the more realistic value of Sum — 1.0 Jy (see Section 2) for source detection we find 
bps — 20-^22 X 10^^ fiK'^SY^ in the Q band. On the other hand, by integrating the T098 model 
counts at 61 GHz (WMAP V band) up to Sum = 0.75 Jy we obtain ^ 1.3 ^ 1.4 x 10^^ 
liK^sr"^, showing a reduced scatter of ~ 0.80 with the Komatsu et al. (2003) best fit value 
(see their Table 1). Moreover, Hinshaw et al. (2003) have found an excess angular power 
spectrum Q = (15.5 ± 1.7) x 10"^ iiK'^sr at small scales, which they explain as produced by 
undetected point sources, in agreement with the findings of Komatsu et al. (2003) . Again, 
by adopting Sum =0.75 and 1.0 Jy for source detection, our current estimates at 41 GHz 
give C( ~ 20 and 24.0 nK'^si, respectively, showing a similar offset as for the point source 
bispectrum. Again, in the V band we obtain Ce 4.5 ^ 5.0 x 10^^ fiK'^sT with Sum = 0.75 
in much better agreement with the Komatsu et al. (2003) estimate of (4. 5 ±4) x 10^'^ fiK'^sr. 
As a result, and in agreement with Komatsu et al. (2003) we actually find an offset between 
their current results and our estimates of bpg. Anyway, Komatsu et al. (2003) results and 
our current estimates - obtained by directly integrating the T098 differential counts - are 
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Fig. 2. — Comparison of tlie reduced CMB bispectrum, in terms of ~ i'^bm, for /„; = 10, thick 
continuous line, with the estimates of bps given by Poisson distributed sources at 30, 44 and 70 
GHz and with Sum = 1 (thick dashed, point-dashed and dotted Hues, respectively). We also plot 
the result with Sum = 0.1 Jy, thin hues. 
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consistent at the la level and, moreover, in the V band the offset (~ 0.8 0.9) is smaller 
than in the Q band. 

Finally, it is important to remind that all these estimates have been computed by 
adopting a Poisson distribution of sources in the sky and a standard CDM cosmology with 
JIm = 1-0. Anyway, our current result are not very much affected by the cosmology: in fact, 
if we use a fiat ACDM model (Qa = 0.7), our estimates only increases of a factor ~ 1.1 1.25 
in the 15-150 GHz frequency interval and by a slightly greater factor {c^ 1.5 ^ 2.3) at higher 
frequencies, while keeping the same parameters as in T098 for the cosmological evolution of 
sources. 



3. THE BISPECTRUM DUE TO CLUSTERED POINT SOURCES 

T098 found that the clustering contribution to CMB fluctuations due extragalactic 
sources is generally small in comparison with the Poisson one. On the other hand, if clustering 
of sources at high rcdshift is very strong it can give rise to a power spectrum stronger than the 
Poisson one (Perrotta ct al. 2003; Scott and White 1999). We remind, however, that current 
models (Toffolatti et al. 1998; Guiderdoni et al. 1998) suggest a broad redshift distribution of 
extragalactic sources contributing to the CMB angular power spectrum, implying a strong 
dilution of the clustering signal (Toffolatti et al. 1999). Anyway, in view of extracting 
the most precise information from current and future CMB data, i.e. pursuing precision 
cosmology, it is also important to assess to which extent the clustering signal can reduce the 
detectability of the primordial CMB bispectrum. 

In a companion paper (Toffolatti et al. 2003) we have addressed the problem of esti- 
mating the CMB angular power spectrum given by undetected clustered point sources by 
adopting a simple phenomenological approach. In that paper we exploit the very few data 
available so far on the angular correlation function of radio and far-IR selected sources to 
estimate the of extragalactic sources at microwave frequencies. We use here the same 
approach and the same simulated maps to estimate the clustering contribution to bps as 
previously done for the Poisson case. We focused on WMAP and Planck LFI channels, 
at which extragalactic radio sources, displaying a "fiat" energy spectrum (Bennett et al. 
2003b), dominate the bright counts. We also give a first estimate of bps at 545 GHz (a 
Planck HFI channel) and we defer to a future paper (Gonzalez-Nuevo et al. 2003) for a most 
comprehensive multifrequency analysis. 

We briefly remind here the adopted procedure. We have carried out 100 simulations in 
2D-flat patches of the sky of 12°.8 x 12°.8 and 25°.6 x 25°.6 deg^ and with a pixel size of ~1.5 
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arcmin for Planck and of ~6 arcmin for WMAP, since maps are created in the HEALPix 
format with nside = 512 (Gorski et ah 1998) in this latter case. A different pixel size can 
highly increase/decrease the pixel-by-pixel variance in presence of a strong clustering signal, 
and the effect has to be taken into account. In this case, at WMAP and LFI frequencies, the 
effect proves negligible in agreement with the previous findings of Toffolatti et al. (1998). 

As discussed in more detail elsewhere (Toffolatti et al. 2003), we first distribute point 
sources at random in the sky with their total number fixed by the integral counts. A^(> S), 
of the T098 model at each given frequency. We then calculate the Fourier transform of the 
density contrast map, i.e. the map defined by 6{x) = ^^^^^^^ , where (A^) is the average 
number of sources per pixel, obtaining a "white noise" power spectrum, P{\<.)poiss =< l^feP 
thus normalized to the total number of sources foreseen by the model. The next step has 
been to modify this normalized spectrum by the Fourier tranform of the angular correlation 
function, w{9), corresponding to the appropriate source population. We adopt here the w{9) 
measured at 4.85 GHz by Loan, Wall and Lahav (1997), since the source population of bright 
sources {S > 50 mJy) which has been sampled at this frequency is representative also of 
bright sources seen at higher frequencies: the underlying parent population appears to be 
the same one (Bennett et al. 2003b) in agreement with the assumption of the T098 model. 
Other estimates of w{9) for radio sources come from samples at lower frequencies (e.g., 1.4 
GHz) and down to lower fluxes where the dominant source populations are different from 
the ones relevant for WMAP and Planck observations. 

With this choice of 'u;(6') the modified density contrast of each pixel in the Fourier space 
is 



where P(k)ci is the Fourier transform of the chosen w{6). In this way the source density is 
modified by the adopted correlation function whereas the normalization to the total effec- 
tive number of sources predicted by the model remains fixed by (N). By the 2D Fourier 
antitransform we get again the map in the real space (the sky patch) . As discussed in detail 
by Toffolatti et al. (2003), the procedure is "safe" since the total number of sources and 
the number counts remain unchanged whereas the recovered w{6) matches very well the 
input one. We want to stress here that this is a phenomenological approach, since we are 
currently interested in reproducing the observational angular correlation functions of the 
relevant source populations, if they are available. 




(12) 
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Finally, from the simulated maps, in terms of AT/T, we then obtain the Fourier transform 

a{k) = J ^^e-^'^'^dh (13) 

and the calculation of the bispectrum in the 2D fiat-sky approximation is given by 

< a(ki)a(k2)a(k3) >= i27rf5^^\kr + k^ + kg) 6(ki, k^, kg) (14) 

where a(k) is the Fourier transform of the temperature map, 5*^^-' is the 2D Dirac delta 
function and 6(ki, k2, ks) is the reduced bispectrum. 

This formula is given by the collapse of the Gaunt integral, Eq.(8), to the Dirac delta 
function in the fiat approximation. We then estimate the reduced bispectrum by means of 
the following estimator, similar to the one used by Santos ct al. (2002), 

Khh = Tf- — -Re[a(ki)a(k2)a(k3)] (15) 

where Re means the real part (which guarantees that the estimator is real). Finally, we 

average over the combination of modes with multipoles li, I2, 13, which satisfies the condition 
ki + k2 + k3 = 0, taking into account Eq.(14). The plotted values refer always to the reduced 
equilateral bispectrum, i.e. for h = I2 = h- 

In Figure 3 we plot the reduced equilateral bispectrum, hps-, calculated at 30 GHz and 
with source detection limits of Sum =1,0.1 Jy. As discussed before, we calculate the average 
values over 100 simulations and the corresponding ilc confidence intervals. To avoid the 
overlapping of too many lines we decided to only plot the confidence intervals corresponding 
to the WMAP resolution given that the ones calculated at the Pianclc resolution are very 
close to the former ones. From the comparison of the top panels (correlated distributions) 
with the bottom ones (Poisson distribution), and with the direct estimates of hps, by inte- 
grating the differential counts (Eq.(lO)), we can extract the following conclusions: a) from 
the simulations of Poisson distributed sources we can estimate hps values with a very high 
precision. The scatter around the straight horizontal line is always much smaller than the 
±1(T limit, b) if we perform simulations with a correlated distribution of sources, e.g. by the 
w{9) of Loan, Wall and Lahav (1997), the value of hps remains practically unchanged. This 
also happens if we choose another angular correlation function among realistic ones which 
could be suitable for the underlying source population. As previously discussed, this is due 
to the very broad redshift distribution of sources contributing to the bright counts in this 
frequency range. On the other hand, the standard deviation results higher for the correlated 
distribution if compared to the Poisson case and, in particular, at low I. In Figures 4 and 5 
we plot our results on hps at WMAP and Pianclc LFI frequencies. It is apparent that there 
are very small variations in comparison with the values of Table 1. 
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Fig. 3. — Reduced equilateral bispectrum, bps, at 30 GHz for two different source detection limits. 
In each panel, the straight horizontal line represents the value of Table 1., obtained by directly 
integrating Eq.(lO) and, thus, using a Poisson distribution of sources in the sky. For comparison, we 
also plot the bps values obtained by averaging over 100 simulations of clustered point sources, upper 
panels, and of Poisson distributed sources, lower panels, at the Planck and WMAP resolutions (see 
text). In the case of WMAP we actually performed the simulations at 33 GHz. The it la limits 
refer to the set of simulations at the WMAP resolution limit. 
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Fig. 4. — hps due to a clustered distribution of point sources at four WMAP frequencies: 23, 41, 
61 and 94 GHz. As in Figure 3, the straight horizontal line represent the h^g values of Table 1. The 
dashed line represents the average 6ps over 100 simulations and the dotted lines represent the ztlcr 
limits. 
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On the other hand, to give a prehminary estimate of the bps at Planck HFI frequencies, 
we have also calculated the reduced equilateral bispectrum of a non Poisson distribution 
of sources at 545 GHz. At this frequency the source populations dominating the number 
counts are dust-enshrouded elliptical galaxies and spheroids at substantial redshift (Toffolatti 
et al. 1998; Guiderdoni et al. 1998; Granato et al. 2001) and, consequently, the clustering 
term should not be negligible (Magliocchetti et al. 2001; Perrotta et al. 2003), given that 
the dilution of the signal is less effective than that of sources showing a broad redshift 
distribution. The results plotted in Figure 6 are obtained, as before, by distributing sources 
in the sky using a Poisson distribution which is then modified by the angular correlation 
functions appropiate for each source populations. As at lower frequencies, the value of bps 
estimated by Eq.(15), applied to the simulated maps and with a Poisson distribution of 
sources, still gives a very good approximation of the value directly obtained by Eq.(lO). 

As for the correlated distributions, for radio selected "flat" -spectrum sources, which are 
giving a quite small but not negligible contribution to the bispectrum at these frequencies, 
we have still applied the w{9) of Loan, Wall and Lahav (1997), under the assumption that 
the underlying parent population is the same (fiat-spectrum QSOs, blazars) up to these 
frequencies. As for sources whose emission is dominated by cold dust, the lack of direct 
data at 545 GHz forced us to rely on SCUBA as well as on optical and near-lR surveys (see 
Perrotta et al. (2003)). By distinguishing the relevant populations in starburts/spirals at 
low redshift and ellipticals and spheroids at intermediate to high redshift, we then applied 
to them the w{6) of Tegmark et al. (2002) and of Maghocchetti et al. (2001), respectively. 
This is a very preliminary estimate and, naturally, it will be improved in the next future. 
Anyway, given that the most relevant quantity is the total number of sources below the 
detection threshold, the choice of a different - but realistic - angular correlation function 
does not modify very much the estimated bps'. Figure 6 clearly show that the three different 
evolution models for sources are giving basically the same estimated bps- On the other hand, 
it is also apparent that at frequencies > 300 GHz source clustering can greatly enhance the 
value of bps and, particularly, at multipoles i <100. This last result is again in agreement 
with the results of Scott and White (1999); Magliocchetti et al. (2001); Toffolatti et al. (2003) 
on the angular power spectrum. 
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Fig. 5. — bps due to a clustered distribution of point sources at Planck LFI frequencies: 30, 44, 70 
and 100 GHz. The meaning of the lines is the same as in Figure 4. 



-20- 



4. DISCUSSION AND CONCLUSIONS 

The primordial bispectrum is a telling quantity of the non-Gaussianity level of the 
CMB and its detection should be of great importance for theories of the early Universe. 
On the other hand, undetected extragalactic sources, which are not Gaussian distributed 
in the sky are surely giving rise to a bispectrum detectable by current as well as future 
CMB experiments. In fact, Komatsu et al. (2003) already claimed the first detection of 
the bispectrum due to undetected point sourecs at 41 GHz by the analysis of the first 
year WMAP data. Therefore, the possible detection of the primordial CMB bispectrum is 
severely hampered by the presence of the foreground emission of point sources which has to 
be carefully evaluated. 

Assuming that the sources are Poisson distributed in the sky their reduced bispectrum 
is easy to calculate by integrating the differential counts of the relevant source populations 
at a given frequency. In this paper, we have carried out a detailed analysis of the bispectrum 
due to extragalactic point sources in all Planck and WMAP frequency channels. We used the 
Toffolatti et al. (1998) evolution model for radio and far-IR selected sources as a template to 
estimate the bps foreseen in the case of a Poisson distribution of point sources and by using 
two different cosmologies (see Section 2). As proved by KSOl, the point source bispectrum 
can be separated from the primordial one taking into account the difference in shape. Our 
current est imates proves that bps can be detected in all WMAP channels if the detection 
limit for sources is Sum >1 Jy ^ albeit at the 1 ^ 2a level at 94 GHz - and also in all 
Planck channels if Sum >0.1 Jy. Table 1 and Fig. 1, 2 show in detail how the values of bps 
of Poisson distributed sources can affect the detection of the primordial bispectrum, bm, in 
the whole frequency interval around the GMB intensity peak and at Planck LFI frequencies, 
respectively. 

In the case of Poisson distributed sources, our main results can be summarized as follows: 

a) the bispectrum detected by Komatsu et al. (2003) in the WMAP 1-year sky maps 
(Q and V bands) is compatible at the la level with the predictions on bps calculated by the 
number counts of T098. The best fit value of bp,^ measured by Komatsu et al. (2003) in the 
Q band can be explained by undetected extragalactic sources according to the T098 model 
predictions multiplied by ~ 0.70, if we integrate the number counts up to Sum =0.75 Jy like 
in Komatsu et al. (2003). On the other hand, in the V band we find a smaller mismatch 
(~ 0.8 -7- 0.9) with current WMAP measurements. Both uncertainties in the measured bps 
and errors in T098 predictions could explain the detected offsets. In particular, the different 
correction factors which are currently found for the Q and V bands, if confirmed by future 
analyses, could be indicative that the adopted spectral index distribution of sources has to 
be partially corrected. Current surveys at 15 and 31 GHz (Waldram et al. 2003; Bennett et 
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Fig. 6. — bps at 545 GHz and applying the conservative source detection limit Sum = 1 Jy (see 
text). As in the previous Figures, the straight horizontal line represents the value of Table 1, 
obtained by Eq.(lO): please note that in the two panels the plotted value is the same one, i.e., 
bps = 5 X 10~^^, in spite of the great scale difference. The top panel shows the comparison of 
the average values of bps obtained by performing 100 simulations with the number counts of three 
different evolution models (see label on the Figure). We adopted different clustered distribution 
of sources for each one of the three different models of sources counts (sec text). The lower panel 
displays the comparison of the analytical estimate of bps with the result of 100 simulations obtained 
by a simple Poisson distribution of sources (dashed line, T098 model counts). In the two panels, 
the dotted lines represent the ±la limits obtained with the T098 source counts. 
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al. 2003b) are showing a greater fraction of "steep" -spectrum radio sources at intermediate 
to bright fluxes than that foreseen by the original T098 model. These new results can help 
in reducing the observed offsets between current measurements of bps and model predictions 
while keeping a good fit to the CBl, VSA and WMAP differential source counts. Another, 
albeit small, source of uncertainty is the underlying cosmology: correction factors of 10-20% 
in the estimated bps values are easily introduced in this frequency range by changing Qm and 

b) The primordial CMB bispcctrum, bm, appears higher than bps only at multipoles close 
to £ ~ 200 where it reaches its peak value (sec Figures 1 and 2). The frequency window inside 
which the peak value of bm is detectable sets larger by increasing and by reducing the 
source detection limit, Sum, i.e. it strongly depends on the efficiency of the source detection 
algorithm. If /„; = 100 and with Sum — 0.1 Jy, the window spans from ~ 30 to ~ 500 
GHz. On the other hand, if fni — 10 and Sum — 2 Jy, the window shrinks to 100 < u < 350 
GHz. These results clearly show that bps does not hamper the possible detection by Planck 
of very low values of the coupling parameter (/„; < 10), comparable to the theoretical limit 
achievable by the mission. 

c) If the detection of bps is effectively achieved (e.g., by WMAP data) the simple com- 
parison with the bps levels estimated by the integration of the differential counts of sources, 
dN/dS, allow, in principle, to test the adopted cvohition model for sources. Anyway, as 
already noted by KSOl, this can be applied only to the bright end of the number counts, at 
least for frequencies below 150^200 GHz, given the 'ffat' slope of the counts - close to or 
(3 = 2.5, i.e. 'euclidean' - which implies that faint fluxes are not relevant in Eq.(lO). As also 
discussed in Section 2., the adopted cosmology does not affect all the above conclusions. 

For estimating to which extent the clustering of sources can affect the value of bps we 
have carried out flat 2D simulations in sky patches with pixel sizes of 1.5 (Planck) and 
6 (WMAP) arcmin. In this analysis we simulated sources with a particular choice of the 
angular correlation function, but we neglected higher order moments of the distribution. 
This simplified approach could affect, in principle, the conclusion that source clustering does 
not modify the bispectrum at frequencies < 100 GHz (see Figures 3, 4 and 5). However, 
as discussed in Section 2 and 3, the results are much more determined by the number of 
undetected sources in the sky (the differential counts of bright sources) than by the choice 
of the correlation function or by the inclusion of higher order moments in the simulations. 
We checked that if we apply a different w{6) to the same number counts, i.e. to the same 
evolution model for counts, the estimated bps keeps practically unchanged. 

On the other hand, at frequencies u >300 GHz, the number of sources in the sky at a 
given flux limit is greatly enhanced by the rising energy spectra due to the emission of the 
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cold dust components. Correspondingly, if sources do cluster down to very low flux limits, 

the clustering term can dominate the Poisson one and bps values result greatly enhanced 
(Figure 6). Therefore, it is surely of great astrophysical interest to study the clustering 
properties as well as the cosmological evolution of extragalactic sources in this frequency 
range for obtaining a better assessment of their contribution to the CMB angular power 
spectrum and bispectrum. 

We wish to thank an anonymous referee for his comments and suggestions. We acknow- 
ledge partial financial support form the Spanish MCYT under projects ESP2001-4542-PE 
and ESP2002-04141-C03-01. FA and LT also wish to thank the EC Research Training Net- 
work contract n. HPRN-CT-2000-00124 for partial flnancial support. JGN acknowledges 
a FPU fellowship of the Spanish Ministry of Education (MEC). We thank B. Guiderdoni 
and G.L. Granato for providing us with their model counts of sources at 545 GHz. We also 
acknowledge very fruitful discussions with J.L. Sanz, E. Martinez-Gonzalez, L. Cayon, R.B. 
Barreiro and P. Vielva. 



-24- 



REFERENCES 

Barreiro, R. B., Hobson, M. P., Lasenby, A. N., Banday, A. J., Gorski K. M., and Hinshaw, 
G., 2000, MNRAS, 318, 475 

Barreiro, R. B., Martinez- Gonzalez, E., and Sanz J. L., 2001, MNRAS, 322, 411 

Bennett, C. L., et al., 2003, ApJ, 583, 1 

Bennett, C. L., et al., 2003, ApJ, in press (astro-ph/0302208) 
Bromley, B. C, and Tegmark, M., 1999, ApJ, 524, L79 

Cayon, L., Sanz, J. L., Martinez-Gonzalez, E., Banday A. J., Argiieso, F., Gallegos, J. E., 
Gorski K. M., and Hinshaw G., 2001, MNRAS, 326, 1243 

Cayon, L., Martinez- Gonzalez, E., Argiieso, P., Banday, A. J., and Gorski, K. M., 2003, 
MNRAS, in press. 

Coble, K., et al., 2003, ApJ, submitted (astro-ph/0301599) 

Palk, T., Rangarajan, R., and Srednicki, M. 1993, ApJ, 403, LI 

Perreira, P. G., Magueijo, J., and Gorski, K.M., 1998, ApJ, 503, LI 

Gangui, A., Lucchin, P., Matarrese, S., and MoUerach, S., 1994, ApJ, 430, 447 

Gangui, A., and Martin, J., 2000, Phys. Rev. D, 62, 103004 

Gonzalez- Nuevo, J., Toffolatti, L., and Argiieso, P., 2003, in preparation. 

Gorski, K.M., Hivon, E., and Wandelt, B.D., 1998, in Evolution of Large Scale Structure: 
from Recombination to Garching, p. 

Granato, G. L., Silva, L., Monaco, P., Panuzzo, P., Salucci, P., De Zotti, G., and Danese, 
L., 2001, MNRAS, 324, 757 

Guiderdoni, B., Hivon, E., Bouchet, P., and Maffei, B., 1998, MNRAS, 295, 877 

Hinshaw, G., et al., 2003, ApJ, submitted (astro-ph/0302217) 

Kofman, L., Blumenthal, G. R., Hodges, H., and, Primack J. R., 1991, in Large-Scale struc- 
tures and peculiar Motions in the Universe, edited by Latham, D. W.& da Costa, L. 
N., ASP Conference Series, Vol.15, p.339 



-25- 



Komatsu, E., and Spergel, D. N., 2001, Phys. Rev. D, 63, 063002 

Komatsu, E., Wandelt, B. D., Spergel, D. N., Banday, A. J., and Gorski, K. M., 2002, ApJ, 
566, 199 

Komatsu, E., Kogut, A., Nolta, M. R., Bennett, C. L., Halpern, M., Hinshaw, G., Jarosik, N., 
Limon, M., Meyer, S. S., Page, L., Spergel, D. N., Tucker, G. S., Verde, L., WoUack, 
E., and Wright, E. L., the WMAP team, 2003, ApJ, in press (astro-ph/0302223) 

Loan, A. J., Wall, J. V. & Lahav, O., 1997, MNRAS, 286, 994 

Luo, X. 1994, ApJ, 427, L71 

Magliocchetti, M., Moscardini, L., Panuzzo, P., Granato, G.L., De Zotti, G., and Danese, 
L., 2001, MNRAS, 325, 1553 

Mandolesi, N., et al., 1998, proposal submitted to ESA for the "Planck" Low Frequency 
Instrument (LFI) 

Mason, B.S., et al., 2003, ApJ, in press, (astro-ph/0205384) 

Pando, J., Valls-Gabaud, D., and Fang, L. Z., 1998, Phys. Rev. Lett., 81, 4568 

Park, C. G., Park, C., and Ratra, B., 2002, ApJ, 568, 9 

Perrotta, F., Magliocchetti, M., Baccigalupi, C., Bartelmann, M., De Zotti, G., Granato, 
G.L., Silva, L., and Danese, L., 2002, MNRAS, 338, 623 

Polenta, G., et al., 2002, ApJ, 572, L27. 

Puget, J.-L., et al., 1998, proposal submitted to ESA for the "Planck" High Frequency 
Instrument (HFI) 

Pyne, T., and Carroll, S.M., 1996, Phys. Rev. D, 53, 2920 

Salopek, D. S., and Bond, J. R., Phys. Rev. D, 42, 3936 

Santos, M. G., et al., 2002, Phys. Rev. Lett., 88, 241602 

Scott, D., and White, M., 1999, Astronomy & Astrophysics, 346, 1 

Schmalzing, J., and Gorski, K. M., 1997, MNRAS, 297, 355 

Sokasian, A., Gawiser, E., and Smoot, G. F., 2001, ApJ, 562, 88 



-26- 



Spergel, D. N., Komatsu, E., Kogut, A., Bennett, C. L., Halpern, M., Hinshaw, G., Jarosik, 
N., Limon, M., Meyer, S. S., Page, L., Nolta, M. R., Tucker, G. S., Verde, L., Wollack, 
E., and Wright, E. L., the WMAP team, 2003, ApJ, submitted (astro-ph/0302209) 

Taylor, A. G., Grainge, K., Jones, M. E., Pooley, G. G., Saunders, R. D. E., and Waldram, 
E. M., 2001, MNRAS, 327, LI 

Taylor, A. C., et al., 2003, MNRAS, in press, (astro-ph/0205381) 

Toffolatti, L., Argiieso Gomez, P., de Zotti, G., Mazzei, P., Pranceschini, A., Danese, L. and 
Burigana, G., 1998, MNRAS, 297, 117 

Toffolatti, L., Argiieso, P., De Zotti, G., and Burigana, C., 1999, in ASP Gonf. Series Vol. 
181, Microwave Foregrounds, A. de Oliveira-Gosta and M. Tegmark eds., p. 153 

Toffolatti, L., Gonzalez-Nuevo, J., and Argiieso, P., 2003, in preparation. 

Vielva, P., Martinez- Gonzalez E., Gallegos, J. E., Toffolatti, L., and Sanz, J.L., 2003, MN- 
RAS, in press (astro-pli/02 12578) 

Waldram, E., Pooley, G.G., Grainge, K.J.B., Jones, M.E., Saunders, R.D.E., Scott, P.P., 
and Taylor, A., 2003, MNRAS, in press (astro-pli/0304275) 

Wang, L., and Kamionkowski, M. 2000, Phys. Rev. D, 61, 63504 

Wu, J. H. P., Balbi, A., Borrill, J., Perreira, P. G., Hanany, S., Jaffe, A. H., Lee, A. T., 
Oh, S., Rabii, B., Richards, P. L., Smoot, G., Stompor, R., and Winant, G.D., 2001, 
ApJS, 132, 1. 



This preprint was prepared with the AAS IM^jX macros v5.0. 



